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RADIATION  PATTERN  OF  RAYLEIGH  WAVES  FROM  A  FAULT  OF 
ARBITRARY  DIP  AND  DIRECTION  OF  MOTION  IN  A 
HOMOGENEOUS  MEDIUM 

By  N.  A.  Haskell 
ABSTRACT 

Expressions  for  the  displacements  in  the  body  waves  radiated  in  an  unbounded,  homogeneous 
elastic  medium  by  dipolar  point  sources  of  arbitrary  orientation  may  be  readily  derived  in 
Cartesian  coordinates  from  formulae  given  by  Love.  The  free-surface  boundary  conditions 
are,  however,  most  conveniently  expressed  in  terms  of  Sezawa’s  cylindrical  wave  functions. 
The  necessary  transformation  between  the  two  representations  is  provided  by  the  Sommerfeld 
integral  and  others  that  may  be  derived  from  it  by  differentiations  with  respect  to  the  radial 
and  axial  (vertical)  coordinates.  By  this  means  the  total  radiation  field  (direct  plus  surface 
reflected)  is  expressed  in  terms  of  integrals  of  cylindrical  wave  functions.  The  Rayleigh  wave 
component  may  then  be  separated  out  by  calculating  the  residue  at  the  Rayleigh  pole  of  the 
integrand.  The  azimuthal  dependence  of  the  Rayleigh  wave  displacements  appears  as  the  sum 
of  three  terms,  one  independent  of  the  azimuth  angle,  <f>,  another  depending  upon  sin  <f>  and 
cos  <f>,  and  a  third  depending  upon  sin  2 <f>  and  cos  2#.  The  coefficients  of  these  terms  are  func¬ 
tions  of  the  direction  cosines  of  the  normal  to  the  fault  plane  and  the  direction  of  the  relative 
displacement  vector  in  the  fault  plane.  Equations  are  presented  ;or  sources  of  both  single  and 
double  couple  types.  The  effect  of  fault  propagation  with  finite  velocity  over  a  finite  distance 
may  be  included  by  multiplying  these  expressions  by  the  finite  source  factor  previously  derived 
by  Ben-Menahem. 

Polar  plots  of  the  amplitude  and  initial  phase  are  presented  for  single  and  double-couple 
representations  of  a  number  of  different  types  of  faults.  It  is  noted  that  for  one  certain  orien¬ 
tation  a  shallow  double-couple  source  generates  no  Rayleigh  waves. 

Introduction 

Ben-Menahem  (1961)  has  recently  treated  the  radiation  of  seismic  surface 
waves  from  moving  single-couple  sources  representing  dip-slip  and  strike-slip- 
faults.  The  case  of  arbitrary  direction  of  motion  in  the  fault  plane  can  be  treated  by 
resolving  the  displacement  into  dip-slip  and  strike-slip  components  and  superim¬ 
posing  the  solutions  for  these  two  special  cases,  but  it  is  also  possible  to  derive 
general  expressions,  allowing  the  normal  to  the  fault  plane  and  the  direction  of 
motion  in  the  fault  plane  to  be  arbitrary  orthogonal  vectors  from  the  beginning, 
which  will  exhibit  the  dependence  of  the  radiation  pattern  on  these  vectors  some¬ 
what  more  explicitly.  We  shall  begin  with  a  representation  of  the  elastic  wave 
radiation  field  of  point  sources  in  terms  of  spherical  wave  functions  in  Cartesian 
coordinates  and  then,  by  means  of  the  Sommerfeld  integral  and  its  derivatives, 
transform  this  into  a  representation  in  terms  of  cylindrical  wave  functions.  The 
free-surface  boundary  conditions  are  then  imposed  and  the  Rayleigh  wave  com¬ 
ponents  are  separated  by  taking  the  residues  at  the  Rayleigh  pole. 

Cartesian  Representation  of  Dipolar  Point  Sources 

In  tensor  notation  the  elastic  wave  displacement  vector  with  Cartesian  com¬ 
ponents  Ui ,  due  to  a  point  force  with  components  F)(t)  may  be  written  (Love, 
1944)  in  the  form 
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pRlfi 

4:7t pv*i  =  i£~3(3fi7i  ^0*)  /  tfFj(t  —  t')  dt'  +  7* y j Fj( i  —  R/a)/aR 

J  r  I*  (1) 

—  (yi  yj  —  5ij)Fj(t  —  R/0}/fi*R 

The  convention  of  summation  over  repeated  dummy  indices  is  to  be  understood 
and  the  following  notation  is  used : 
p  =  density 

a  =  compressional  wave  velocity 
0  =  shear  wave  velocity 

7 i  =■  direction  cosine  of  line  from  source  to  field  point 
=  (xi  “  x/)/R 

Xi  —  Cartesian  coordinate  of  field  point 
Xi  —  Cartesian  coordinate  of  source 
R  —  radial  distance  of  field  point  from  source 
dij  =  0,  i  t*  j 

1,  %  j 

Following  Keilis^Borok  (1950),  equation  (1)  may  be  written  in  a  more  symmetrical 
form  by  introducing  the  second  integral  of  the  force,  J (/),  i.e.,  F(t)  =  The 

integral  may  then  be  evaluated  and  equation  (1)  becomes 

47 rpUi  -  (37*7 i  —  ~~  R/<x)/R 3  +  (37*7;  ~  —  R/ct)./aR 2 

+  7,7 -  R/a)/o?R  -[(37,7y  -  ~  R/P)/R3  (2) 

+  (37,- yj  -  in)  J  fit  -  R/(t)/0l  +  (7,7,-  -  hi)  J  fit  -  R/fi)/f?R ] 

To  obtain  the  displacements  due  to  dipolar  point  sources,  we  need  the  first  deriv¬ 
atives,  Uitk  =  dui/dxk  .  From  the  definitions  we  have  dR/dXk  =  7* ,  dyi/dXk  = 
(&ik  —  7*7/0 /R,  and  dJj(t  —  R/a)/dXk  —  —y *«//(<  —  R/a)/a.  Using  these  rela¬ 
tions  we  find 


4:TTpUitK  —  J  j(t  —  R/a)[S(yi8jk  +  7  fine  +  —  157t7;7/J/^4 

+  J/(t  —  R/a)[Z(yiSjk  +  yfiik  +  y»iu)  —  157*7  3yk\/<xR? 

+  J  /'(t  —  R/a)[yi8jk  +  yjSik  +  yk&u  —  Qyz'YjyK]/aR2 
-  J/"(t  -  R/a)ytyjyk/aR 

(3) 

“  J —  R/ff)  [3(7*5 jk  +  yfiik  +  7 k&ij)  —  157t7;7/fc]/-^4 
—  J  /(t  —  R/ /?)  [3  (yudjk  +  yfiik  +  7/ Ai)  —  157*7j7J/^3 
~  J j,r  ( i  —  ^/^)[7*^;A  +  y j^ik  4*  27it6*y  —  67*7 SYkI/I^R2 
+  J/"(t  —  R/fi)[yty3yk  —  7 k$ii\/&zR 
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We  now  assume  that  Jj  =  —  (fi/of)  exp  (icd) ,  where  fj  is  a  unit  vector  in  the 
direction  of  Fj ,,  so  that  Fj  —  Jf  —  fj  exp  {fat)  for  a  couple  of  unit  moment. 
Omitting  the  common  time  factor,  equation  (3)  becomes 

4nrpuiik  =  fjR~l  exp  (  —  io)R/a)  [(3/Rz)  {y$3k  -4  7y$#  *+*  7 k8ij  —  5yz'yjyk) 

+  (Sfa/aR?)  {y:i8  3k  +  yfiik  +  7 kh?  "  57<y  jjlk) 

+  («2/ ctR) (Oy^Yjy.k  —  Vidjk  —  YjSik  ~  7 khij)  +  yj  (4) 

+  //R  1  exp  {—faR/p)  [(3/ R*)(yidjk  +  yjdjk  +  7^*7  ~~  57*7*7*) 

+  ( Sfa/pR2)  (7 $jk  +  yj&ik  +  yjihj  67*7*7*) 

+  (co2 /p2Ry(6y <r/y *  “  7Ak-  VA*  —  27^)  +  (fa*/Pz)  (7*7/7*  -  7*$o')l 

Let  rc*  be  a  unit  vector  normal  to  the  fault  plane.  We  shall  consider  only  shear 
faults,  i.e.,  no  component  of  displacement  normal  to  the  fault  plane,  so  that  rn  and 
fi  are  orthogonal  vectors.  The  displacements  due  to  a  single-couple  source  whose 
orientation  is  specified  by  the  vectors  and  fi  are  then 

Uj  =  —Ui'kUh  (.5) 

From  equation  (4)  and  the  orthogonality  condition,  n^fi  =  0,  the  displacements 
reduce  to 

4 mni  -  Fa[niCf{^/Rz)  +  (3 fa/atf)  -  (J/ aR)}. 

+  f>Cn{(3/R3)  +  (3 w/atf)  -  (<S/c?R)} 

-  y  iCfC„{  (15/ R3)  +  (15m/ aR1)  -  (6a f/dR)  -  («//«*)}]; 

(6) 

-  FAniCAWR?)  +  (3 w/0)  -  («*/ tfR)} 

+  /iC,{(3 /R3)  +  (3 ia,/PRr)  -  (2 J/tfR)  - 

-  %CfCnt(15/R3)  +  (15m/ ISR2)  -  (6 u/fR}  -  (m /!?)]) 


where 

Fa  =  RT1  exp  (—icoR/ a) 

F$  =  RT1  exp  (—ibiR/R) 

(7) 

Cf  =  yji  ■ 

C„  =  Y*W* 
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Transformation  to  Cylindrical  Coordinates  and  Wave  Functions 

Equation  (6)  is  now  to  be  transformed  into  a  cylindrical  coordinate  system 
(r,.  zr  <t>)  with  origin  at  the  source. 

r  -  Hxt  -  xt')2  +  (x2  ^  x2')*\112 

z  =  Xz  —  Xz 

cos  <j>  =  (xi  —  xx)  /r 

In  terms  of  these  coordinates  and  the  spherical  distance,  R,  the  direction  cosines 
%i  are 


71  =  (r/jR)  cos  <t> 

72  =  ( r/R )  sin  <£  (8) 

73  =  z/R 

and  the  cylindrical  components  of  the  displacement  (ur ,  uz ,  u$)  are 

ur  =  (R/r)(yiUi  +  y  2t^) 

Uz  =  uz  (9) 

u*  =  (R/r)(yiU2  -  y2Ui) 

The  scalar  products  Cf  and  Cn  become 

Cf  =  R~l(rfi  cos  <j>  +  rft  sin  <f>  +  fzz) 

(10) 

Cn  ~  R  (mi  cos  <j>  +  m2  sin  +  n3z) 

Using  the  orthogonality  of  and  m ,  their  product  may  be  written 
C/Cn  ==  (t  /2R  )t(/i%  /2W2)  cos  24  +  (fi7i2  +  fffli)  sin  2 <f>\ 

+  (^/^2)[(/i%  +  fan  1)  cos  cj>  +  (/27i3  +  /3w2)  sin  <t>]  (11) 

% 

+  fM(z/R)2~  (r/R)2/2\ 

It  will  be  convenient  to  collect  terms  according  to  their  dependence  on  the  azimuth 
angle,  <j>,  by  writing 

a 


UT  =  Wr(0>  +  Url)  -f  Urm 


(12) 
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and  similar  expressions  for  uz  and  % ,  where  the  terms  with  superscript  (0)  in¬ 
clude  those  independent  of  0,.  the  (1)  term  are  those  depending  on  sin  <j>  and  cos  0 
and  the  (2)  term  are  those  depen  ding,  on  sin  2 0  and  cos  20.  When  this  separation  is 
carried  out,  it  is  found  that  the  coefficients  of  the  various  angular  functions  can 
be  expressed  in  terms  of  the  derivatives  up  to  the  third  order  of  the  functions 
Fa  and  Fp  with  respect  to  r  and  z,  for  example, 

dFa/dr  =  —(rFa/'BHJT1 4  io>/«) 


d-Fa/drdz  =  {rzFa/R2)[{ 3/fi2)  4  (3 iu/aR)  -  (»*/«)! 
etc.  The  algebra  is  straightforward  but  tedious  and  we  give  only  the  final  results: 
4t rpurm  =  (fm/2nSd%/drdz*  4  (co/ a)1  (dFa/dr)  -  3 d%/drdz2],  (13) 

4irp«,(1)  =  cos  -4  j/n1).(d3Fa/dr’dz  —  d3FJ:drdz)  —  j\n3(u>/fi)3dFp/dz} 

(14) 

4  sin  0[(/2w3  4'  fm)(d3Fa/dr2dz  —  d3Fp/dr2dz)  —  f2n3(a/P)2dF.p/dz] 
4?rpMr<2)  =  (l/2)[(/j«x  —  /2n2)  cos  20  4  (fm+fcni)  sin  2<t>][2dsF„/dr* 


(15) 

4  d3Fa/drdz2  4  (c i/a'fdFJdr  -  2 d3Fp/dr3  -  dzFp/drdz*  —  2 (w/&)2dFp/dr] 


4jrpw/0)  =  (/3M3/2)  \3d3Fa/dz3  4  (o>/ct)2dFa/dz  —  3d3Fp/dz3  —  3(w/i S)2dFp/dz] 
4jrpMi(1)  =  cos  0f(/i«3  4  ftfh)  (d3Fa/drdz2  —  d3Fp/drdz2)  —  f3ni(u/p)2dFp/dr] 
4  sin  0[(/2n3  4  /3M2)  (d3Fa/drdz2  —  d3Fp/0rdz 2)  —  /3«2(o)//3)  2dffy/dr] 
47rpMs(2)  =  (1/2)  [(/i«i  —  /2W2)  cos  20  4  (/i«2  4  /*Wi)  sin  20][2d3/i,a/dr2da 
4  d3Fa/ilz3  4  (u/afdFJdz  -  2d3Fp/drdz  -  d%/dz3  -  (w/0)2dFp/dz] 
4t p%(0)  =  (1/2)  (/i%  -  /2W1)  (a/fifdFp/dr 

iirpu^  =  cos  0[(/2%  4  /3W2)  (d2Fa/rdrdz  —  d‘F p/vdrdz)  —  f2n3(ai/^)2dFp/dz\ 
—  sin  0[(/i«3  4  /3W1)  (d2Fa/rdrdz  —  d2Fp/rdrdz)  —  fin3(o>/$)2dFp/dz] 
47rptt/2i  =  (l/2)f(/jfi2  4  /2wi)  cos  20  —  (/pix  —  /2W2)  sin  20][— 2d3Fa/dr3 
—  2d3 Fa/ drdz2  —  2  (  a/ a)  2dF a/ dr  4  2 d3Fp/ dr3  4  2d3Fp/  drdz2  4  (w/jd)  5/ dr] 


(16) 

(17) 


(18) 


(19) 

(20) 


(21) 


The  above  expressions  for  the  displacements,  although  expressed  in  cylindrical 
coordinates,  still  involve  the  spherical  wave  functions  F„  and  Fg .  In  order  to 
complete  the  transformation  into  a  cylindrical  wave  representation  we  use  the 
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Sommerfeld  integral  (Ewing,  Jardetsky,  and  Press,  1957  pp  13  and  14)  for  Fa 
and  Ffi 


where 


Fa  =  R-1  exp  (- iuR/a ).  =  /  J*(kr)e~tM{k/va)  .1 

Jq 

f+V^ft2  ( o)/a )2 

«  I  , _ 

( 4*i\/ (w/a)2  —  fc2 


(22) 

/c  >  w/qj 
k  <  w/a 


Calculating  the  indicated  derivatives  and  using  the  recursion  relations  for  the 
Bessel  functions,  the  various  displacement  components  for  z  >  0  are  as  follows : 


47rpwr(0)  =  (fznz/2)J  Ji(kr)[3vp  e~vpz  — •  (3pa  +  oy  /  ava)e~VaZ]k2  dk 

4irptt,<0)  —  {fznz/2)  j  Jo(kr)[{2(oi/a)2  —  3 k2}e~VaZ  +  3 k2e~V(iZ}k  dk 

4t rpw^<0>  =  (1/2)  (f2m  -  /xn2)(w/iS)2  [  Ji(kr)e~VfiZ(k2/vfi)  dk 

Jo 

4:Tpura)  =  cos  ^  j~(/i«3  +  jjf  Jo(kr)  (e~y“z  —  e~Vf>z)k3  dk 

Ji(kr)(e~’“‘  —  e~y^)k2  dkj  ~h  finzioi/P)2  J  Jt(kr)e~"f‘k  d/cj 
+  sin  </>  [(*„  +/3W2)  jjf  Jo(kr)(e~raZ  —  e~v?z)k3  dk 
~r  l{  Jl(kr^e~"°‘  ~  e-^)fc2dfoj  +/2w3(co//3)2jf  J0(kr)e-yn  rf/cj 
4»pMia)  =  cos  4>  +  /3W1)  jf  Ji(kr)(vpe~'’^  —■  va e~"aZ)k2  dk 

+  hniiWP)2  J [  Ji(kr)e-^(k2/pfi)  dfcj 
+  sin  <j>  j^(/2W3  +  f?,rh)  jf  Ji(kr)(vpe~r$z  —  vae~v*z)ki  dk 
■  fznaiw/P)2  Jilkrie^itf/ve)  dfcj 


(23) 

(24) 

(25) 


(26) 


(27) 


+  . 


J 
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iirpu^l)  =  cos  <l>  £(/2«3  +  /3n^)r_1jf  Ji(*rK*”  ~  dk 

+  Ms(a>//3ffo  Jo{hr)e^n ;  <&J  -  sin  <j>  jy,*,  fM.Jr'1 

•  f~  Ji(kr)(e~’“  -  e~vntfdk  +  /msCw/tffj [  Jo(kr)e~^k  dfcj. 

4fl7>wr(2)  =  (l/2)[(/ini  —  /2R2)  cos  2^5>  +  (/i^2  +  /2%i)  sin  20] 

*  f  Ji(kr){e~v<xZ(k*/va).  —  die 

Jo 

—  (2/r)  f  Jz(kr)  \va~1d~y<x!!  —  j^rV***}^* dk 
4 

47r/>w2<2)  =  (l/2)[(/i7ii  —  72^2)  cos  20  +  (/i^2  +  /2&1)  sin  20] 

*  f  JtlkrXf"  -  <Tv**)kzdk 

4 


(28) 


(29) 


(30) 


4?rpw/2)  —  (l/2)[(fifl*  +/2W1)  cos  20  —  (/ltti  — /2W2)  sin  20] 

•f(2/r>  f*  /tflfcr)C*V’*  -  dk  (31) 

+  («/£)*  f  JiOcr)e~^(k2M  dk 
Jo 

The  corresponding  expressions  for  z  <  0  may  be  obtained  by  changing  the  sign 
of  va  and  v$  and  then  reversing  the  sign  of  all  terms. 

These  expressions  have  now  been  converted  into  the  form  of  integrals  of  the 
solutions  of  the  elastic  wave  equations  in  cylindrical  coordinates  in  the  form  given 
by  Sezawa  (1931)  .  Sezawa’s  cylindrical  wave  functions  may  be  written  in  the 
form 


4 rpurln)(k)  =  cos  nj,lFinc(e)kJn-i(kr)  -  (n/r)Jn(kr){Finc(z)  +  F3nc(z)}] 
+.si i#i“(i)WMW  -  (n/r)Jn(kr)iF1n,(z)  +  F2n’(z)}] 
4?r ptir(n)(k)  =  cos  n^-Fi\z)Jn{kr)  +  sin  n^-F/*(z)  Jn(kr) 

4 W#)(*0  =  sin  vtWrWJ^kr)  -  {n/r)Jn(kr)\Fr{z)  +  ^3nc<0)}] 
-  cos  nbWrmJ^ikr)  -  (n/r)Jn{kr){Ft‘(z)  +  Fzn,{z))) 


(32) 


(33) 


(34) 
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The  functions  Fin(z )  and  F2n(z)  can  be  written  in  the  form 

Fxn(z)  =  Z2n'(z )  -  V(s)  (35) 

F2B(s)  =  k% n(z)  -  Zi(s)  (36) 

where  Z\n  and  Z2  are  solutions  of 

Zi  -  v«%n  =  0  (3.7). 

Z2n  -  ve%n  »  0  (38) 

The  function  .F$n{z)  satisfies  the  same  differential  equation  (38)  as  Z2{z). 

We  now  take  the  origin  of  coordinates  at  the  source  with  the  +z  axis  pointing 
vertically  downward  into  the  medium.  Let  the  depth  of  the  source  be  h,  so  that 
the  free  surface  is  at  z  —  —  h.  In  the  region  above  the  source,  —h  <  z  <  0,  both 
upward  and  downward  travelling  waves  exist  and  the  appropriate  solutions  of 
equations  (37)  and  (38)  are  of  the  form 

Zic  -  Acnle~PaS-  +  A°n 26™ 

Zr  =  Asnle~VaZ  +  Asn2eVaZ 

zr  =  +  Bl  2ev*z 

(30) 

Zf8  -  Bsnle~^  +  B8n2ev** 

Fznc  -  Ccn,e^z  +  Ccn2ev&z 
Fr  -  C*nle~v&z  +  Cl*’* 

The  values  of  the  coefficients  of  the  terms  in  eVaZ  and  e^z  may  be  determined  by 
identifying  terms  in  the  integrated  solution 


=  f  urM(k)  dk 

Jo 


with  the  corresponding  terms  in  the  counterparts  of  equations  (23)  through  (31) 
for  2  <  0.  The  results  of  this  are: 

A  02  =  kf3n:i[2(o}/af  —  3fc2]/2v„ 

BI2  =  —  3/cf3«3/2 

C02  =  (fori!  —  /jn2)  </k/2l?ve 

•A  12  —  k*  (fans  -f-  fan r)  A  Is  =  1c  (fitis  +  hnf) 

Bn  =  finsVp  +  hnjc/ye  Bn  =  /2m3^  +  (41) 

=  firk(o>/B?  Cl 2  =  f2n3(a>/0f 

Alt  =  (/2w2  -  /j r?i)fc3/2j'a  ii22  =  —  (/i«2  4-/2%)&3/2ra 

B22  =  —  (fith  A-  hnf)k/2 
Cn  =  —  (/i«2  +  fon^uk^fvf. 


Bn  =  (/2w2  -  fini)k/2 
C22  =  (/2M2  ~  finx)(ak/2 /J2 
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The  coefficients  A°ni ,  etc.  are  then  determined  by  the  vanishing,  of  the  stress 
components  Tzz ,  Tzr ,  Tg<f>  at  the  free  surface.  These  are 

Tzz  =  XA  +  2 pdus/dz 

Ttr  =  p(duz/dr  +  dUr/dz)  (42) 

Tz<f>  =  fi^du^/'dz  duz/rd<j>) 

where  A  =  dur/dr  +  ur/r  +  dujrd<j>  +  duz/dz  and  X  and  p  are  the  Lam4  elastic 
constants.  Setting  these  stress  components  equal  to  zero  at  z  —  —ft  leads  to  the 
following  set  of  equations 

(1  -  7 XAV*  +  AUT+)  +  yp^-BUe^  +  K^)  =  0 
-  (7  -  1  )W>h  +  BbT*)  -  (yvJk*)(-AW“h  +  Al*f*h)  =  0  (43) 

-Cl  1  +  Cl,e~^h  =  0 

where  7  =  2(/3/c/co)2,  and  a  similar  set  for  the  coefficients  with  superscript  s.  Solving 
for  Ant ,  etc.  we  have 

Acnl  =  [-^(A)4c„2e-2-“A  +  27(7  -  D 

Bh  =  [{27(7  ~  -  ff(*)B»2e-2^]//(A)  (44) 

<Tnl  =  e~2^hCCn2 

where 

/(*)  -  (y-  l)2-  t2v/ 

?(*)  -  (t  ~  l)2  +  yW*2 


Radiation  Pattern  of  Rayleigh  Waves 

The  part  of  the  displacement  field  due  to  Rayleigh  waves  is  the  contribution  to 
the  integral  of  equation  (40)  that  arises  from  the  singularity  of  the  integrand  at 
the  real  root  of  /(ft).  By  a  transformation  of  the  path  of  integration  in  the  complex  k 
plane  (Ewing,  Jardetsky,  and  Press  1957,  pp  132-135)  this  reduces  to  the  negative 
of  the  residue  at  the  pole,  k  =  k,  where  /(/c)  =  0.  This  may  be  obtained  by  mul¬ 
tiplying  the  right  hand  sides  of  equations  (32),  (33),  and  (34)  by  —  7r i,  replacing 
the  Bessel  functions  by  the  corresponding  Hankel  functions  of  the  second  kind, 
Hu2)(ht),  dropping  those  terms  of  FVi  F2 ,  and  Fz  that  do  not  contain  f(k)  in  the 
denominator,  and  replacing  f(k)  by/(/c)  in  the  denominator  of  the  remaining  terms. 
At  this  point  we  shall  also  drop  the  “near-field”  terms  in  (n/r)Hn{2)  (kr)7  since 
they  are  of  order  r~d/2  and  thus  make  a  negligible  contribution  to  the  distant  radia- 
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tioii  field  in  comparison  with  the  leading  terms  which  are  of  order  r  1/2.  The  re¬ 
maining  significant  terms  are  then 

4 pUrn)  =  —  ulln-i (kt) [Gi”c (z)  cos  ntj)  +  G^8(z)  sin  n^] 

4 pU2in)  ==■  —  @»(2,)  W cos  n<j>  +  £2™  (2)  sin  n<£]  (45) 

Vtt)  -  0 

where  the  functions  £?(z)  are  derived  from  the  functions  F(z)  by  applying  the 
operations  described  above.  Using  equations  (35),  (36),  (39),  (41),  (44),  together 
with  ./(*)  =  Q,  explicit  expressions  for  the  values  of  these  functions  at  the  free 
surface  (z  =  —h)  are  as  follows: 

Gi°c  =  {^yu^n3/f(K)]{{(y  -  l)/y}fh’>  -  {1  -  2c#V3«*} «-*'*]  (46) 

GiQe  =  yvaGic/(y  -  1)  (47) 

G\C  =  [2k  (y  -  1 )//(<)] 

(48) 

•[-{/3W1  +  /i«3(l  -  2/7)6  +  (y  —  1)  (fins  +  ftni)e  kVa/y] 

Giu  =  [2  k2  (7  -  1  )  //  (  *)] 

(49) 

•[-{/3W2  +  /s«»(l  —  2/7)  e  +  (7  —  l)(/2«3  +  /3w2)e  *”a/7] 


(?2lc  =  yvaGi"/  (7  -  1)  (50) 

(?2ls  =  yvaGx‘ / (y  -  1)  (51) 

<?i2c  =  [k7^(M  -  /2«2)//(k)][((7  -  DA}*"*"  -  «~**1  (52) 

ft2*  =  [w(/i%  +  /2«i)//(K>]t{(7  -  l)/i)e-^  -  e-k’°]  (53) 

(?22c  =  yvaG?°/(y  -  1)  (54) 

G ?  =  yvaGi‘/(y  -  1)  (55) 

Replacing  the  Hankel  functions  by  their  asymptotic  approximations 

Hn2y(ia')  — »  (2/tkt)112  exp  (— i/cr  +  (inn/  2)  +  zV/4)  (56) 

the  surface  displacements  reduce  to 
ur  =  4(k,  r,  ft)e-w/4[3*/3n3{(l  -  2cRl/Zot)D  -  (y  -  l)/y] 


+  {2k(7  —  l)/yvp]  cos0[— {/3%  +/i?i3(l  —  2/7)}  +  (7  —  1)  (/ift3  -H/3rei)Z)/7] 

(57) 

+  {2k(7  -  D/yvfi]  sin <t>[— [fin2  +f2n3(l  —  2/y}}  H-  (7  —  1)  (Ms-\- M2)D/y] 

4-  i[{  (7  —  1)/ 7}  —  D][(/i«i  -  M2)  cos  2<t>  +  (M2  +  Mi)  sin  2<f>)  ] 

uz  =  iyvauT/n(y  —  1)  (58) 


RADIATION  PATTERN  OF  RAYLEIGH  WAVES 

A(k,  r,  h)  “  {K27^/4p/'(«:)}{2/7r/cr}1/2  exp  (—ur  —  hvp) 
D  =  exp{-h(va  -  vp)\ 
cR  =  u/tc.  “  phase  velocity  of  Rayleigh,  waves 


In  interpreting  the  phase  angles  of  the  displacements  in  terms  of  the  direction 
of  motion  on  the  fault,  we  note  that  the  normal  vector,  ny  may  be  drawn  on  either 


Fig.  1.  Convention  for  relation  between  /  and  n  vectors  and  fault  displacement. 

TABLE  1 

Values  of  /i  Corresponding  to  Various  Types  of  Faults 


Type 

h 

h 

Dip  slip 

Strike  slip 

Dip  slip 

Strike 

slip 

Dip  slip 

Strike 

slip 

Normal,  right-lateral 

0 

—1 

COS  5 

0 

sin  8 

0 

Normal,  left-lateral 

0 

1 

COS  8 

0 

sin  8 

0 

Reverse,  right -lateral 

0 

-1 

■  **-cos  8 

0 

—sin  5 

•  0 

Reverse,  left-lateral 

0 

1 

—  COS  8 

0 

—sin  8 

0 

side  of  the  fault.  The  sign  convention  adopted  in  equation  (5)  then  requires  that 
the  vector  /  be  interpreted  as  the  direction  of  F  on  the  same  side  of  the  fault  as  that 
on  which  n  is  drawn.  Since  in  the  immediate  neighborhood  of  the  source  the  direc¬ 
tion  of  F  is  the  same  as  that  of  the  displacement,  the  relationships  between  /, 
n,  and  the  directions  of  motion  on  the  fault  plane  are  those  shown  in  figure  1.  If 
we  choose  the  x\  axis  to  be  along  the  strike  of  the  fault  and  the  xo  axis  in  the  direc¬ 
tion  of  dip,  with  the  directions  chosen  as  shown  in  figure  2,  and  let  n  be  drawn  on 
the  hanging  wall  side  of  the  fault,  the  components  of  n  are  nx  —  0,  n2  —  sin  5, 
ns  —  —cos  8  and  the  limiting  values  of  the  components  of  f  for  various  types  of 
faults  are  those  given  in  table  1. 

Double  Couple  Source 

It  must  still  be  considered  an  open  question  whether  the  majority  of  earthquakes 
are  best  represented  by  a  single-couple  source  model  or  by  a  double-couple  with 
zero  resultant  moment.  To  obtain  the  surface  wave  radiation  pattern  for  the  latter 
case  we  note  from  figure  1  that  interchanging  the  vectors  n  and  /  produces  a  couple 
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at  right  angles  to  the  original  couple  with  moment  of  opposite  sign.  The  appro¬ 
priate  solution  for  the  double-couple  source  model  may  then  be  obtained  from  the 
single-couple  solution  by  interchanging  corresponding  components  of  n  and  f 
and  adding  the  results  to  the  original  solution.  The  result  is 

ur  =  2 A(k,t,  ~  2c„2/3a2)Z»  -  (7  -  DM 

4  2 (vtt/k)  (D  -  1){  (/in,.  4  fnh)  cos-  <t>  4  ( fiit, 3  4  /*n»>  siii  <t> }  (61) 


-  i{D  -  (7  -  1)/7}.{(/i%  -  fnn2)  cos  2 <j>  4  (/iw2  4/«%>  sin2M  J 


Fig.  2.  Orientation  of  axes  with  respect  to  fault. 

Where  we  have  made  use  o//(k)  =  0  in  deriving  the  coefficient  of  the  terms  in  sin  <f> 
and  cos  <j>.  The  relation  between  uz  and  ur  continues,  of  course,  to  be  that  given 
by  equation  (58) . 


Effect  of  Fault  Propagation 

Ben-Menahem  (1961)  has  treated  the  effect  on  the  surface  wave  radiation  pat¬ 
tern  of  horizontal  propagation  of  the  fault  fracture  at  finite  velocity  over  a  finite 
distance.  Although  the  case  of  horizontal  propagation  of  the  fracture  is  probably 
the  most  common  one,  his  method  can  be  readily  extended  to  cover  the  general 
case  of  an  arbitrary  direction  of  propagation.  Let  £  be  the  distance  from  the  point 
of  initiation  to  the  instantaneous  position  of  the  leading  edge  of  the  fracture  meas¬ 
ured  in  the  direction  of  fracture  propagation,  and  let  v  be  the  velocity  of  fracture 
propagation.  Let  n ,  zq  ,  0O  be  the  cylindrical  coordinates  of  the  point  of  observa¬ 
tion  with  respect  to  the  point  of  initiation  of  the  fracture,  and  r,  z}  <j>  be  the  corre¬ 
sponding  coordinates  with  respect  to  the  instantaneous  position  of  the  leading  edge 
of  the  fracture.  Then,  following  Ben-Menahem  we  may  write 


—  fb 

u(ro ,  zo ,  <f>o)  =  F1  /  u(r,  z,  <j>)  exp  {iw(t  —  %/v)  d£\ 
vo 


(62) 


) 

i 


u 


1 
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where  b  is  the  length  of  the  fault,  u(r0  zq ,,  0O)  stands  for  any  one  of  the  displace¬ 
ment  components,  and  u(r9  zr  <f>)  is  any  one  of  the  previously  derived  expressions 
for  a  stationary  point  source.  Let  (fa  ,  fa ,  fa, ,  be  the  Cartesian  components  of  a 
unit  vector  in  the  direction  of  fault  propagation.  Physically  l  must  be  normal  to 
n  so  that  Itfii  ==  0,  but  otherwise  its  direction  is  unrestricted.  Then  to  terms  of  the 
first  order  in  £  we  have 

r  ==  ?’o  —  £(A  cos  0  +  fa  sin  <f>)  -f  0(£2) 

£  =  Zo  “  &  (63) 

0  —  <j>6  +  (£/r)  (fa  sin  (j)  ~  fa  cos ,<j>)  -f*  0.(£2) 

Now  consider  the  azimuth  dependent  factor  in  u(r ,  z ,  4>)  .  This  contains  terms  of 
the  form  cos  n<f>  and  sin  mj>  and  we  have 

cos  Ticf)  =  cos  n<f>o  +  (0  —  <f) 0)  d  cos  n0o/d0 0  +  •  •  • 

=  cos  m&o  —  (£/r)  (4  sin  <f>  —  4  cos  0)n  sin  n0o  +  0[(£/r)2] 

with  a  similar  expression  for  sin  n<f> .  Thus  if  we  consider  the  radiation  pattern  at 
distances  such  that  r  ^S>  6,  so  that  terms  of  order  £/r  are  negligible,  we  may  treat 
the  azimuthal  terms  in  the  integrand  as  constants  in  carrying  out  the  integration 
in  equation  (62).  Similarly  for  the  radial  amplitude  factor,  r_1/2,  we  have 

r“1/2  =  nT1/2  +  (r  —  ro)  dnll2/dro  +  •  •  • 

—  ro  1/2[1  +  (£/2ro)>(2i  cos  0  +  ?2  sin  0)  +  0{(£/ro)2}] 

and  this  factor  may  also  be  treated  as  a  constant  to  the  same  order  of  approxima¬ 
tion.  Therefore  the  only  factors  that  need  to  be  considered  in  evaluating  the  integral 
in  equation  (62)  are  the  exponential  radial  and  depth  factors.  Noting  that  at  the 
free  surface  z  =  —h  and  z0  =  —  h0 ,  so  that  h  =  k o  +  the  only  factors  that 
that  need  be  retained  under  the  integral  are  those  of  the  form 

/  exp  (—inr  —  hvp  —  ifc/v)  d£  =  exp  (—  uro  ~  hovp) 

Jo 

(64) 

•  /  exp  [iK%(l i  cos  (f> o  +  fa  sin  0o )  —  £fa  vp  —  zco£/V|  d£ 

Jo 

A  similar  expression  with  a  replacing  /3  also  occurs  through  the  quantity  D  that 
appears  in  the  coefficients  of  the  azimuthal  terms  in  equations  (57)  and  (61). 
Apart  from  a  small  change  in  the  relative  amplitudes  of  the  azimuthal  terms  due 
to  this  quantity,  the  main  effect  of  source  propagation  is  to  multiply  the  expressions 
for  stationary  sources  by  the  factor 
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where 


a  —  u(li  cos  <j>  +  k  sin  <j>  —  cR/v )  —  kvp,  (66) 

The  subscript  indicating  the  origin  at  the  initial  point  of  fracture  has  now  been 
dropped.  In  the  case  of  horizontal  propagation  of  the  fault,  k  =  0,  and  if  we  set 

X  =  (jtb/2y{(cR/v)  —  lx  cos  (f>  l2  sin  <j>\  .(67). 

equation  (65)  reduces  to  the  form  given  by  Ben-Menahem 

P(k,  </>)  =  (sin  X/X)  exp  (WZ).  (68) 


Initial  Phase 

The  phase  angles  of  the  displacements  ur  and  uz  obtained  by  setting  r  =  0  in  the 
expressions  developed  above  are  the  initial  phases  in  the  sense  of  the  term  employed 
by  Aki  (1960  a,  b,  c)  and  Brune  (1961,  1962).  They  do  not,  of  course,  represent 
the  actual  phases  of  the  ground  motion  near  the  epicenter,  since  they  do  not  in¬ 
clude  the  near-field  terms  that  make  a  negligible  contribution  at  large  distances, 
but  are  the  dominant  terms  at  short  distances.  The  initial  phase  of  the  far-field 
terms  does,  however,  have  physical  significance  since  it  is  the  phase  angle  that 
one  obtains  by  applying  the  phase  compensation  methods  of  Aki  or  Brune  to  ob¬ 
served  data  obtained  at  large  distances. 

Since  the  expressions  given  above  refer  to  a  source  with  sinusoidal  time  de¬ 
pendence,  an  additional  initial  phase  factor  is  necessary  if  these  expressions  are  to 
be  applied  to  a  particular  Fourier  component  of  a  nonsinusoidal  wave  form.  If  the 
time  dependence  of  the  force  couple  representing  the  source  is  F(t),  with  a  Fourier 
transform 


/+00 

F{t)e~iutdt 

•oo 

the  correction  that  should  be  added  to  the  initial  phase  is  arg  [/(«)].  In  particular, 
if  F(t)  is  a  unit  step  function,  /( a>)  =  1/fco,  and  arg  [jf(o>)]  =  —tv/2  should  be 
added  to  the  initial  phase  of  all  components. 

In  the  case  of  Rayleigh  waves  on  a  homogeneous  medium  the  quantity  is 
real  and  positive.  From  the  definition  of  f(k)  and  making  use  of  the  fact  that  k 
is  the  value  of  k  such  that  /(*)  —  0,  it  may  be  shown  that 

/'(*)  =  ~27[72(1  -  M  +  (1/y)  ~  2]A(7  -  l)2  (69) 

Since  y  =  2 and  cR  <  /?  for  all  homogeneous  media,  y  >  2.  Also  a  = 
/z/(A  +  2fi)r  which  is  necessarily  <§.  Hence  72(1  —  a )  >  2  and  the  term  in 
brackets  is  a  positive  real  quantity.  Therefore  /'(«)  is  a  negative  real  quantity  and 
the  factor  A  (k,  r,  h)  in  equation  (57)  has  a  phase  angle  of  rr. 

For  direct  comparison  of  initial  phase  angles  computed  from  equations  (57), 
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(58),  and  (61)  with  the  initial  phase  in  the  sense  of  Brune,  two  additional  correc¬ 
tions  are  needed.  Brune  takes  the  positive  sense  of  the  2-axis  as  directed  upward 
instead  of  downward  in  accordance  with  the  convention  that  we  have  used  in  the 
present  paper.  This  corresponds  to  a  change  of  =fc.x  in  the  phase  angle  assigned  to 
uz .  He  also  makes  a  correction  of  — x/4  to  remove  the  phase  advance  at  the  source. 
In  the  present  theory  this  phase  advance  is  represented  by  the  factor  exp  (irf/ 4) 
that  occurs  in  the  asymptotic  expression  for  the  Hankel  function,  equation  (56) 
and  must  be  removed  for  consistency  with  Brunei  definition. 

In  summary,  if  we  let  &•(<£)  be  the  phase  angle  of  the  factor  in  brackets  in  equa¬ 
tions,  (57)  or  (61)  that  expresses  the  azimuthal  dependence  of  the  initial  phase  of 
ur ,  the  phase  angle  corrections  that  must  be  added  to  produce  Brunei  initial 
phase,.  j  are  as  follows: 


TABLE  2 


Summary  of  Illustrative  Cases 


8- 

M2 

m 

h 

h 

f 9 

Single-couple 

Double-couple 

90 

1.0 

0.0 

-1.0 

0.0 

0.0 

Fig.  3 

Identical  with  fig.  3 

85 

0.9962 

-0.0872 

-1.0 

0.0 

0.0 

Fig.  4 

— 

80 

0.9848 

-0.1736 

-1.0 

0.0 

0.0 

Fig.  5 

— 

70 

0.9397 

-0.3420 

-1.0 

0.0 

0.0 

Fig.  6 

— 

45 

0.7071 

-0.7071 

:  -1.0 

0.0 

0.0 

Fig.  7 

Identical  with  fig.  3  with 
scale  reduced  by  0.7071 

0 

0.0 

-1.0 

^1.0 

0.0 

0.0 

Fig.  8 

Amplitude  —  0 

90 

1.0 

0.0 

-0.7071 

0.0 

0.7071 

Fig.  11 

Identical  with  fig.  3  with 
scale  reduced  by  0.7071 

90 

1.0 

0.0 

0.0 

0.0 

1.0 

Fig.  12 

Amplitude  —  0 

45 

0.7071 

-0.7071 

:  -0.7071 

0.5 

0.5 

Fig.  13 

fig.  14 

45 

0.7071 

-0.7071 

0 

0.7071 

0.7071 

Fig.  15 

fig.  16 

Phase  correction  to 
be  added  to  ^r(<£) 

1)  exp  ( —  7rt/4)  factor  in  eqs.  (57)  and  (61)  —  x/4 

2)  Negative  j sign  of  /'(«)  in  A(k,  0,  h) 

3)  Step  source  function  —  r/2 

4)  Conversion  to  vertical  component  {uz  ~  iur)  +ir/2 

5)  Change  to  opposite  sign  convention  of  2-axis  rfcx 

6)  Phase  advance  at  source  —  7r/4 


Net  correction  —  x/2 

so  that  ypB(4>)  —  ^v(<£)  V2 

Representative  Examples 

Radiation  patterns  have  been  computed  for  a  number  of  illustrative  cases  and 
are  summarized  in  table  2,  The  strike  of  the  fault  plane  is  taken  as  the  reference 
direction  (<£  =  0)  so  that  n\  =  0  in  all  cases.  The  direction  of  dip  is  down  to  the 
right  and  all  figures  refer  to  the  case  of  zero  focal  depth  (Z>  =  1.0). 


SINGLE  COUPLE.  RT  LAT.,  STRIKE  SLIP 

Fig.  4.  Right-lateral,  strike-slip  motion  on  a  fault  dipping  85°.  Single-couple  source  mo 


PARAMETER  ;  b/X 

Fig.  9.  Amplitude  of  fault  propagation  factor  for  horizontal  propagation  at  v  =  /3.  X  ^  b» 
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PARAMETER  =  b/X 

Fig.  10.  Amplitude  of  fault  propagation  factor  for  horizontal  propagation  at  v  =  ,0.  X  <  6. 

The  effect  of  fault  dip  on  the  radiation  pattern  for  the  single-couple  source  model 
of  a  right-hand,  strike-slip  fault  is  illustrated  in  figures  3  through  8,  which  show 
polar  plots  of  relative  amplitude  and  initial  phase  (Brune  convention  with  step 
function  source)  as  functions  of  azimuth.  It  will  be  noted  that  the  quadrantal 
symmetry  with  two  nodal  axes  that  exists  in  the  case  of  a  vertical  fault  evolves 
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nithei  rapidly  into  a  two  lobcd  pattern  with  a  single  nodal  axis  as  the  dip  dec  reases. 
These  figures  lefer  to  the  ease  of  a  stationary  source.  To  visualize  the  effect  of 
fault  propagation  the  amplitudes  must  be  multiplied  by  the  Ren-Mcnahcin  propa- 
tion  factoi,  which  is  shown  as  a  function  of  azimuth  in  figures  1)  and  10  for  various 
ratios  of  fault  length,  hf  to  wave-length,  X.  The  fault  propagation  velocity  is  taken 
to  be  equal  to  the  shear  wave  velocity,  so  that  c/v  =  .9194,  and  horizontal  propaga¬ 
tion  in  the  0=0  direction  is  assumed.  For  b/\  <  0.125  the  modification  of  the 
stationary  radiation  pattern  for  amplitudes  is  small  (<10%).  For  b/\  >  .521, 
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SINGLE  COUPLE.  RT  LAT.,  NORMAL,  45#  PLUNGE 

Fra.  11.  Amplitude  and  initial  phase  (Bruno)  radiation  pattern  for  45°  plunging  motion  on  a 
vertical  fault.  Single-couple  source  model. 

additional  nodal  axes  due  to  fault  propagation  are  introduced,  and  as  5/X  becomes 
»1  the  radiation  becomes  increasingly  concentrated  into  two  beams  at  cos  0  ™ 
c/c,  or  0  ■*  =b23.1G°,  For  the  double-couple  source  model  the  Rayleigh  wave  rmlui- 
tion  pattern  is  identical  with  the  single-couple  model  in  the  case  of  strike  slip  on  a 
vertical  fault.  However,  for  dips  less  than  90°  the  double couple  source  continues 
to  give  the  same  fourdobed  pattern,  the  only  difference  being  a  decrease  in  ab¬ 
solute  amplitude  in  proportion  to  the  sine  of  the  dip  angle.  At  zero  dip  the  amplitude 
becomes  zero  (for  zero  focal  depth  only). 

Figures  11  and  12  show  the  patterns  for  the  single  couple  source  model  <  i  a 
vertical  fault  with  the  direction  of  motion  plunging  45°  and  90°  (dip  slip)  respec¬ 
tively,  The  effect  is  the  same  as  changing  the  dip  of  a  strike  slip  fault  except  that 
the  pattern  is  rotated  through  (KF.  Tn  these  cases  also  the  patten  for  the  double 
couple  source  remains  unchanged,  but  the  amplitude  decreases  as  (he  plunge  in¬ 
creases,  and  become^  zero  for  pure  dip-slip  mot  ion. 
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RADIATION  PATTERN  OF  RAYLEIGH  WAVES 

Figure  13  shows  the  single-couple  radiation  pattern  for  45°  dip  angle  and  dis¬ 
placement  vector  making  an  angle  of  45°  with  the  horizontal  in  the  fault  plane. 
Figure  14  shows  the  corresponding  patterns  for  the  double-couple  source.  Although 
the  amplitude  patterns  appear  similar,  the  double-couple  case  actually  lias  a  pair 
of  minor  lobes,  too  small  to  be  shown  in  the  figure,  between  the  two  major  lobes. 
The  differences  in  initial  phase  distribution  would  distinguish  between  these  two 
cases. 

The  corresponding  comparison  between  single  and  double-couple  sources  for  dip- 
slip  motion  on  a  fault  of  45°  dip  is  shown  in  figures  15  and  16.  In  this  case  also 


DOUBLE  COUPLE  f,  *0,  fr  16  fs  =  .7071 

Fig.  16.  Double-couple  source  model.  Orientation  parameters  as  in  figure  15. 


the  amplitude  patterns  are  similar,  but  the  initial  phase  patterns  would  distinguish 
between  the  two  source  models. 
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